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Abstract 

Atmospheric dynamics span a range of time-scales. The projection 
of measured data to a slow manifold, A4, removes fast gravity waves 
from the initial state for numerical simulations of the atmosphere. 
We explore further the slow manifold for a simple atmospheric model 
introduced by Lorenz and anticipate that our results will relevant to 
the vastly more detailed dynamics of atmospheres and oceans. 

Within the dynamics of the Lorenz model, we make clear the rela- 
tion between a slow manifold A4 and the "slowest invariant manifold" 
(sim), which was constructed by Lorenz in order to avoid the diver- 
gence of approximation schemes for A4. These manifolds are shown 
to be identical to within exponentially small terms, and so the SIM in 
fact shares the asymptotic nature of A4 . 

We also investigate the issue of balancing initial data in order to 
remove gravity waves. This is a question of how to compute an "initial- 
ized" point on M. whose subsequent evolution matches that from the 
measured initial data that in general lie off M.. We propose a choice 
based on the intuitive idea that the initialization procedure should 
not significantly alter the forecast. Numerical results demonstrate the 
utility of our initialization scheme. 

The normal form for Lorenz' atmospheric model shows clearly how 
to separate the dynamics of the different atmospheric waves. How- 
ever, its construction demonstrates that any initialization procedure 
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must eventually alter the forecast — the time-scale of the divergence be- 
tween the initialized and the uninitialized solutions is inevitable and 
is inversely proportional to the square of initial level of gravity-wave 
activity. 
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1 Introduction 

Wave motion in the atmosphere has a wide range of periods: large-scale 
motions are dominated by quasi-geostrophic Rossby waves, which have a 
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period of many days; faster inertial-gravity waves, with a period of up to a 
few hours, are frequently negligible in the large (Gill [12, p582]). A significant 
problem for numerical weather prediction is that an initial measured state 
of the atmosphere contains noise which causes unphysical large-amplitude 
gravity waves to arise in numerical solutions (Gill [12, p242]). Consequently, 
one is forced to take a very small time step in the numerical integration 
(e.g. Houghton [16, pl69]). This is despite the generally recognized feature 
that much of the troposphere and lower stratosphere is quasi-geostrophic and 
so dominated by Rossby wave activity (Lorenz [22]). 

To overcome this problem of unphysical gravity-wave activity, initial data 
must be "balanced", that is, adjusted so that subsequent gravity-wave activ- 
ity is negligible (see e.g. Houghton [16, pl68], or Gill [12, p245]). A straight- 
forward balancing procedure is to project the initial data onto the normal 
modes of the system, which depends on the details of the numerical proce- 
dure used for forecasting (Errico [11]); the balancing is achieved by setting 
the gravity wave modes to zero. Such a procedure proves unsuccessful be- 
cause gravity waves immediately appear in the numerical simulations — the 
nonlinear aspects of atmospheric dynamics are large enough to render a linear 
procedure ineffective. Various more sophisticated initialization schemes exist 
(see Daley [10], and references therein), many of which can be interpreted as 
projecting the initial data onto a "slow manifold" (Leith [19]; Lorenz [21]), 
Ai. On such a slow manifold the gravity wave variables are given as functions 
of the Rossby wave amplitudes, and there are no rapid oscillations. 

Lorenz [21, 22] and Lorenz & Krishnamurthy [20] (henceforth referred to 
as L80, L86 and LK87, respectively) have proposed a series of low-dimensional 
models of the atmosphere and find several difficulties with the practical com- 
putation of a slow manifold: that various approximation schemes for Ai are 
divergent; that fast gravity wave oscillations apparently must occur in some 
regions of Ai; and that resonances between gravity waves and Rossby waves 
induce singularities in the manifold. Such complex dynamics of such sim- 
ple models is further explored by Camassa [5] who proves, for example, the 
existence of chaotic dynamics. We concentrate here on one of these mod- 
els (L86), in which, according to the linearized governing equations, neither 
Rossby waves (which we denote by x) nor gravity waves (y) are subject to 
damping. A three-mode Rossby-wave complex is coupled to the two gravity- 
wave modes through nonlinear interactions. The Rossby waves have a much 
longer period than the gravity waves, and in fact the period of the Rossby 
waves in the model becomes infinite as their amplitude tends to zero. 

We emphasize that we analyze only the low- dimensional model of L86. 
The applications of the ideas presented herein to realistic atmospheric equa- 
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tions will involve much more complicating detail. Nonetheless, we antici- 
pate no great difficulty in extending the analysis to more physical dynam- 
ics. The procedures utilized in this paper to construct slow manifolds and 
to appropriately initialize data have already been generalized successfully 
from low- dimensional toys to high-dimensional physical problems (Mercer 
& Roberts [23, 24]; Roberts [28]; and Watt et al [34]). The corresponding 
analysis of actual atmospheric models is left to further work. 

The absence of damping in the model of L86 is a source of difficulties that 
do not arise when the waves are slightly damped, as LK87 considered. One 
might expect therefore that since atmospheric waves clearly are damped, the 
less troublesome case would be the most profitable to explore. Further, the 
limit as the damping a — > is singular, and theoretically the nature of M. is 
different in the two cases a = and a / 0. It might seem perverse, therefore, 
to tackle the undamped case. However, a closer examination reveals that the 
mathematical difficulties do not disappear when damping is introduced; they 
are merely obscured. For if a is assumed to be small then inverse powers of a 
occur in certain calculations, and these are a precursor of the problems that 
arise when a = 0. By examining the undamped case (as does Camassa [5]) we 
confront head-on the fundamental difficulties which result from resonances 
between waves of different timescales. 

We view a rational balancing procedure as having two stages. First, the 
functional form of the slow manifold must be determined — usually in the 
form y = h{x). Once this is accomplished a closed set of low-dimensional 
evolution equations for the slow variables follows. This reduced set could 
be integrated forwards in time to make forecasts, although in practice the 
full system is integrated from the balanced initial data. Secondly, the ap- 
propriate initial values for the slow variables on M. must be calculated from 
the full set of initial values. Once the slow variables x are known, the ap- 
propriate fast variables are then given by y = h(x). The first stage has 
received most attention (Baer & Tribbia [3]; Leith [19]), while the second 
seems relatively ignored in the meteorological literature, although the ap- 
propriateness of some different projections, according to the quality of one's 
data measurements, has been discussed by Daley [9, 10]. In Section 2 we pro- 
pose a different criterion than those previously considered for the selection 
of an initial point on the slow manifold: that the subsequent evolution cor- 
responds as closely as possible to that of the full initial-value problem, but 
with gravity waves absent. (We shall make the notion of "corresponding" 
evolution more precise later.) A similar choice is known in physics, where 
fast variables are eliminated (Haake & Lewenstein [15]; van Kampen [30]), 
and seems to have an obvious desirability in numerical weather forecasting: 
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the adjustment of the initial data should not alter the forecast. Many of the 
approximation schemes to compute the slow manifold implicitly assume that 
only the fast variables require adjustment, although Daley illustrates how 
this is not the case for his "optimal" projection schemes. For our proposed 
projection criterion also, it turns out that both fast and slow variables must 
be altered. 

In Section 2 we briefly describe the formal computation of a slow manifold, 
At, for L86, together with the choice of appropriate initial conditions on Ai. 
In Section 3 we discuss the divergence of the power series for Ai , and consider 
Lorenz' computation of an alternative slow manifold, the "slowest invariant 
manifold" (sim). We find that the SIM and Ai differ by terms smaller than 
any power of x as x — > 0, so that the SIM has a divergent power series. 
A normal form transformation is presented in Section 4, where the model 
L86 is written as the slow evolution of five new variables; the five variables 
include the amplitude and phase of the gravity waves. The normal form 
for L86 has dynamical behavior that differs from that of L86, albeit by an 
exponentially small amount. Although these differences are quantitatively 
asymptotically small, we find that some solutions are qualitatively affected 
by them. Furthermore, the normal form shows that the evolution of the slow 
modes cannot be entirely decoupled from the fast modes. Thus there must be 
inevitable discrepancies in the long-term evolution of the slow waves, of the 
order of the square of the fast waves, between the balanced and unbalanced 
simulations. 

Finally, in Section 5 we review the notions of a fuzzy manifold, and the 
implications of the present work for the initialization of numerical weather 
forecasting schemes. 



2 Quasi-geostrophy as a subcentre manifold 

2.1 Computation of a slow manifold for L86 

In the model of L86, the slow "Rossby wave" variables x = (U, V, W) and 
the fast "gravity wave" variables y = (X, Z) evolve according to 

U = -VW + bVZ 
V = UW-bUZ 



W 



-UV 



(1) 



X 



-z 



z 



X + bUV, 
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where the over-dot denotes differentiation with respect to time, b is a coupling 
parameter, and the superscript T denotes the transpose of a row vector. 
Observe that when these equations are linearized about the zero equilibrium 
the Rossby wave variables remain constant, while the gravity wave variables 
oscillate sinusoidally with a period that has been normalized to 2tt. 

When b = the Rossby waves and the gravity waves are uncoupled: the 
gravity waves oscillate sinusoidally with period 2tt; while the Rossby waves, 
given by Jacobian elliptic functions, oscillate with a period inversely propor- 
tional to their initial amplitude. For small-amplitude motions the Rossby 
waves have a much longer period than the gravity waves — this remains su- 
perficially true when coupling is restored (6^0). 

General solutions of (1) have gravity- wave components; we should like 
to adjust a given initial condition for (1), called balancing, so that gravity 
waves do not develop, while maintaining essentially the same evolution of the 
physically significant Rossby waves. That is, we propose the principle that a 
long-term forecast should be unaffected by the initialization. 

A first attempt to balance initial data might be projection onto the 
"geostrophic manifold", given by y = 0, that is, X = Z = 0. This mani- 
fold is certainly free of gravity waves, but is not invariant under (1); if we 
set X = Z = initially, they do not remain zero. A more sophisticated 
approach to balancing is to seek an invariant "quasi-geostrophic manifold" , 
or slow manifold, M, on which solutions of (1) evolve slowly. Instead of 
seeking y = 0, we allow the fast variables to be functions of the slow vari- 
ables, y = h(x) (Leith [19]). The slow manifold is a so-called subcentre 
manifold (Kelley [18]). Little is known about the conditions for the existence 
of subcentre manifolds; indeed resonances plague attempts to construct such 
manifolds (Sijbrand [29]). In this section we proceed with a formal con- 
struction of A4, leaving the issues of its existence and uniqueness until later 
sections. 

The details of the calculation of M. are given by Roberts [26]. By 
substituting the ansatz y = h(x) into (1), and applying the chain rule 
y = dh/dxx, we obtain a quasi-linear partial differential equation for h(x) 
(Carr [6]). Since this PDE cannot be solved exactly we expand h( X^j cLS Sb 
power series in the slow variables x, and find 

X= -bUV + O (\x\ 4 ) , Z = b{U 2 -V 2 )W + 0(\x\ b ) , (2) 

as | a; | — > 0. Therefore on Ai the Rossby waves evolve according to the 
following (slow) amplitude equations obtained by substituting (2) into (1) 



U = -VW 



l-b(U 2 - V 2 )] +£>(|a;| 6 ) 
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w 



uw 
-uv. 



1 - b(U 2 - V' 



O Ice 



(3) 



The procedure we have used to compute h(x), that is, to find X(U,V,W) 
and Z(U, V, W), is mechanical, and has also been described for the linearly 
damped model L80 by Vautard and Legras [32]. It is equivalent to the 
balancing scheme of Baer & Tribbia [3] . 



2.2 Initialization: projection of initial conditions 

We must supplement our computation of the functional form of the slow 
manifold y = h(x) by deriving appropriate initial conditions for x on M.. 
By "appropriate" we mean that the subsequent evolution on M. faithfully 
follows the behavior of the full system (1) from the original initial conditions. 
In general it is not sufficient to adjust only the amplitudes of the gravity 
waves, that is, to map the initial point (a;*, y*) of the full system to the point 
(cc*, h(x*)) on M.: if the full system and the slow-manifold model are to have 
the same future behavior from their respective initial conditions, the initial 
values of the Rossby wave variables also must be adjusted. The discrepancy 
between the initial condition for the slow variables in the full initial-value 
problem and that on the slow manifold is known as "initial slip" (Grad [13]), 
by analogy with the slip allowed for an inviscid fluid at a boundary. We now 
proceed to calculate this "initial slip". 

The appropriate choice of an initial point on M. has a straightforward 
derivation (Roberts [26]) when M, attracts neighboring solutions exponen- 
tially. In that case, the choice may be made so that the solution from the 
original initial condition, and that from the adjusted initial condition ap- 
proach one another exponentially quickly as t — >• oo. However, in the present 
case the slow manifold does not attract neighboring trajectories, but instead 
acts as a centre for their gravity-wave oscillations — a solution initially off M. 
oscillates about M. perpetually. We therefore aim to choose the initial point 
on M. so that its subsequent evolution maintains its relationship with the 
full uninitialized solution for all time. 

The procedure we describe below is correct to leading order in the dis- 
tance, r, of the initial point (x*,y*) from Ai. In Section 4 we shall show 
that in general it is not possible to improve the choice of initialized point. 
For example, we expected to be able to incorporate corrections of C(r 2 ); 
however, we find that this cannot be done and so solutions on and off M. 
necessarily drift apart after a time of O (r~ 2 ). 
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In order to explain our procedure we begin by briefly considering the 
simple projection 



[x 



y*)^(x*,h(x*)), (4) 

used implicitly by Baer & Tribbia [3] . The adjustment of initial data, by the 
displacement (0, h(x*) — y*), lies in the "gravity-wave" space spanned by 
the vectors (0, 0, 0, 1, 0) and (0, 0, 0, 0, 1): only the gravity-wave variables are 
altered. We call the solution obtained by projecting initial data in this way 
BT. However, for the projection we propose here, the two vectors that span 
the appropriate direction for the projection of initial conditions depend on 
x* (Roberts [26]), and in general the Rossby-wave variables are also changed 
by the initialization process. The simple projection (4) is appropriate only 
in the special case when x* = 0, where we may observe that x(t) = for all 
t > 0, according to both the original system (1) and the slow model (3). 

A more informed projection scheme must be based upon the evolution 
near the slow manifold (Roberts [26]). If U(t) is a solution of (1) on Ai, 
and if e(t) = (ex, e 2 , €3, 64, 65) is an infinitesimal displacement from U then e 
evolves, according to (1), as 



e = (C+Nx)e, 



(5) 



where Ce = (0,0,0,-65,64), and M\ is the Frechet derivative of the vector 
of nonlinear terms H(x) = (-VW + bVZ, UW - bUZ, -UV, 0, bUV) eval- 
uated on Ai (see Roberts [26]). For any slow solution on M., we then seek 
those neighboring solutions that evolve with the solution on Ai plus a fast 
oscillation. It is just these neighboring solutions which should be projected 
onto the particular slow solution on Ai in order to maintain the long-term 
forecast. We wish to calculate only the projection space spanned by ei and 
e 2 , and not their magnitudes nor their individual directions. The argument 
is elaborated in detail in Section 4 of the paper by Roberts [26] . The result to 
low-order is that the projection onto Ai should be made along the position 
dependent planes spanned by vectors ex, e 2 , where 



ei 



-bV 
bU 



1 





+ 



x\ 



e 2 






-b(U 2 - V 2 ) 


1 



O \x 



(6) 



We shall call the solution obtained after projecting initial data along these 
vectors R. 

In Figure 1 we show some typical results for solutions BT and R, and 
compare them with solutions of (1) from uninitialized initial data. We have 
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Figure 1: A comparison* of two ways of projecting thc tl ufitial condition 
(x*,y*) onto A4: BT and R. Each graph shows three solutions of (1): 
the solid line from the initial condition (x*,y*); the dotted line from 
projection R; the dashed line from projection BT. Initial conditions are: 
(i) (x*,y*) = {U,V,W,X,Z) = (0.1,0.1,-0.1,0.01,0.01); (ii) (x*,y*) = 
(—0.01,0.1,0.01,0.01,0). In graph (ii) the solid line and the dotted line are 
almost indistinguishable, except for the small-amplitude gravity-wave wig- 
gles present in the solid line, but not in the dashed. Projection R is clearly 
better than BT. Both initialization schemes use just the first few terms in 
the relevant power-series expansions, as given by (2) and (6). 

chosen different initial conditions for the two graphs; these are representative 
of a large number of numerical runs. In both cases, R is clearly closer to the 
uninitialized solution than is BT. In the first, BT and R are both phase- 
shifted with respect to the original solution, but R less so than BT. Further, 
BT goes the "wrong way" round a fixed point of (1) at time t ~ 50, and 
as a result the slow variable W takes the wrong sign. In the second case, R 
and the original solution are indistinguishable on the figure (although they 
slowly diverge at later times); at times t « 75, 150 the original solution has 
gravity-wave wiggles that the projected solutions do not have (because they 
both lie close to the slow manifold). Here, R is strikingly better than BT. 

2.3 Incorporation of forcing 

The centre manifold formalism allows one to compute the effects of a forcing 
of the full system (1) on the dynamics of the slow manifold (or, for that 
matter, the damped, forced system of LK87). This computation is described 
to leading order in the amplitude of the forcing by Cox & Roberts [7] , and 
shows that the forcing not only changes the evolution of the slow "Rossby" 
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wave variables, but also changes the shape and location of a slow manifold 
M.. For example, if we add a constant forcing F = (Fjj, F v , F w , F x , F z ) to 
the right-hand side of (1) then the slow manifold becomes, to leading order 
in \F\, 

X = -bUV -F Z + {U 2 -V 2 )F W -2W(VF V -UF V ) 

+ 0(\x\\\F\\x\ 3 ) (7) 

Z = bW(U 2 -V 2 ) + F x + b(UF v + VF V ) -b 2 F x (U 2 -V 2 ) 

+ 0(|;c| 5 ,|F||:z| 3 ) . (8) 

We have previously shown (Cox & Roberts [7]) that when forcing is present, it 
is important to project initial conditions onto a slow manifold given by (7-8), 
rather than the unperturbed slow manifold (2), in order to eliminate gravity 
waves. 



3 Non-uniqueness of a slow manifold 

3.1 Divergence of series 

Lorenz [22] demonstrates that several approximation schemes aimed at con- 
structing a slow manifold yield divergent power series. Two of these schemes 
(Baer & Tribbia [3]; Vautard & Legras [32]) are equivalent to the construc- 
tion of a subcentre manifold (Roberts [26]) that we have just described. Not 
surprisingly then, the expansions (2) are found to be divergent. 

Such divergence is not necessarily a computational disaster. Techniques 
such as Pade summation or the Shanks transform not only improve the rate 
of convergence of a series, they may also produce a converged "sum" of a di- 
vergent series — the Stieltjes series is a traditional example (see Section 8.3 in 
Bender & Orszag [4]). The use of such convergence-acceleration techniques 
in fluid dynamics has been reviewed by Van Dyke [31]. Here the series expan- 
sion (2) for A4 is in terms of the three slow variables U, V and W. However, 
the acceleration of convergence of such multi-variate expansions is less well 
understood than that for expansions in a single variable. Perhaps the best 
current technique is to use the multi-variate Pade transform proposed by 
Cuyt [8], but our preliminary investigations are so far inconclusive. At the 
very least, low-order computations of the shape of a slow manifold can be 
summed to within an exponentially small error using an "optimal truncation" 
of the series (Bender & Orszag [4, pp.94-100,122-123]). 
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3.2 Construction of the SIM 

Lorenz (L86) describes an alternative construction of a slow manifold, called 
the "slowest invariant manifold" (sim), for which convergence is not a prob- 
lem. The calculation proceeds by generating a family of periodic solutions 
to (1), using a numerical shooting method to determine the correct initial 
value of one of the variables in order to ensure periodicity. This family forms 
an invariant manifold, the SIM, and the central issue is whether such an ob- 
ject is free of significant gravity-wave activity. The construction is heavily 
dependent on the symmetries of (1). 

Lorenz' method of generating the periodic orbits is to consider solutions 
for which both V* = V^O) and X* = X(0) are initially zero. He then takes 
initial values W* > U* > and treats the remaining initial condition, Z*, as a 
parameter to be chosen so that at the first zero-crossing of U, when t = T say, 
X vanishes too. A consequence of the simultaneous vanishing of U(T) and 
X(T), together with the symmetries of (1) (see L86), is that the constructed 
solution is periodic, with period AT. For a fixed ratio m 1 / 2 = U*/W*, Lorenz 
computes a one-parameter family of periodic orbits, as W* is varied. The 
curve Z = Z*(W*) is then a section through the SIM, and the entire SIM may 
in principle be computed numerically by varying the ratio U*/W*. Lorenz 
finds that the SIM contains an infinite number of singularities, which are 
increasingly closely packed as W* — > 0. He concludes that "there is no 
unequivocally slow manifold" because his candidate slow manifold, the SIM, 
contains regions of high gravity-wave activity close to the singularities. 

The central result of this section will be that Lorenz' SIM and the slow 
manifold M. computed in the previous section are exponentially close, a con- 
cept we elaborate upon below. Furthermore, a slow manifold is not unique: 
indeed, there are infinitely many pretenders to that title. Such nonuniqueness 
arises frequently in the construction of low- dimensional models of dynami- 
cal systems (Roberts [27]), and each has the same (divergent) power-series 
expansion. 

Lorenz shows how the singularities in the SIM arise from resonances be- 
tween the slow Rossby waves and the fast gravity waves, by considering the 
uncoupled system (1) with 6 = 0. Then the section through the SIM in 
the (W*, Z*)-plane consists of the horizontal line Z* = 0, on which there 
is no gravity- wave activity, together with the vertical lines W* = K/(k7r), 
for integers k, on which the period of the Rossby waves is a multiple of the 
gravity-wave period. The singularities of the SIM when b ^ come from a 
perturbation of the structure of the uncoupled system. Lorenz noted that 
the singular curve Z = Z*(W*) looked like Z = aexp(l/W*)cot(K/W*). 
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By computing the analytic form of Z*, up to terms of O (b), we now show 
that this expression for Z* is qualitatively correct, and that the singularities 
are indeed exponentially weak as W* — > 0. To do so we expand the variables 
as power series in the coupling parameter b, so that 

oo 

and so on. (LK87 make a similar expansion for the variables, although with 
a different purpose.) We take initial conditions 

oo 

u(o) = u\ v(o) = o, w(o) = w\ x(o) = o, z(o) = z* ~ v z h 

3=0 

where U* and W* are fixed independently of b, and we let the time of the 
first zero-crossing of U(t) be at t = T = T + bT\ + ■ • -. Recall that in order 
to find periodic solutions of (1), our aim is to find T and Z* so that 

U{T) = X(T) = 0. (9) 

At leading order in b the variables satisfy 

Uo = -Wo 
V = U W 

Wo = -U V (10) 

Xq = —Zq 
Zq = Xq. 

The solution (L86) is 

Uo = U*cn(W*t), V = U*sn(W*t), W = W*dn(W*t), 

X = —Zq sin t, Z = Zq cos t, 

where cn, sn, dn are Jacobian elliptic functions (Abramowitz and Stegun, 
Chapter 16) with parameter m = (U* /W*) 2 . 

Expanding (9) in powers of b, we find at leading order that T satisfies 
Uq(T ) = 0. The first zero-crossing of U (t) occurs when t = T = K/W*, 
where 

K = K{m)=i — - 2 (li; 

1 — m sin 9 
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is the complete elliptic integral of the first kind. Thus X Q (T Q ) = —Zq sin(K/W*), 
which vanishes if either Zq = or if W* = Kj (kir), for some integer k. This 
is Lorenz' result for b = 0. The line Zq = is the flat approximation to the 
SIM, and the vertical lines W* = Kj (kir) indicate resonances between the 
Rossby waves and the gravity waves. We take Zq = 0, so that X = Z = 0, 
and now consider the terms of O (b) in the expansion of the SIM. These 
satisfy 

U x = -VqWx-VxWq 
Vi = U W 1 + U 1 W 

Wi = -UoVi-UM (12) 
X x = -Z x 

z x = x l + u Q v Q , 

subject to 

£M0) = K(0) = ^i(O) = X 1 (0) = 0, Zi(0) = z*. 

We see immediately that U x — V\ — W\ = (and indeed, U 2n +i = V 2n+ i = 
W 2n +i = X 2n = Z 2n = for all n). However, the equations for the gravity 
waves yield 

XAt) = -Zlsint-U* 2 /* sm(i - t) cn(WV) Bn(WV) dr , (13) 

Jo 

and Z x (t) = —X x (t) from (12). Now we examine (9) at O (b), which yields 

T 1 U (T ) = 0, X 1 (T ) = 0. 

It follows that T\ = 0, and upon substituting (13) into the second condition, 
we must choose Z\ so that 1 

Z{ sin T = -U* 2 [ T ° sin(T - r) cn{W*r) sn{W*r) dr. (14) 
Jo 

By contour integration (see the appendix for details) we find 

z; = i sech (pf) cot (w) + W ' F(W1 ' (15) 

as shown in Figure 2, where K' = K{1 — m) and where F(W*) has the 

1 This result may also be derived by considering the terms of O (b) as a forcing of the 
unperturbed system, as in Cox & Roberts [7]. 
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Figure 2: A section through the SIM given by the analytic expression (15). 
asymptotic series 

oo 

F(W*) ~ ^dc (2n) (0|l-m)W* 2n 

n=l 

~ mW* 2 + m(m + 4)W* 4 + m{m 2 + 44m + IQ)W* 6 

+ m(m 3 + 408m 2 + 912m + 64)W* 8 + O (w* 10 ) . (16) 

Here dc^ 2n ' ) is the 2n-th derivative of the Jacobian elliptic function dc. 

Observe that through the factor of cot (if /W*) the behavior of Z\ is 
singular for W* near K/(kir), but there is an exponentially small factor, 
sech(K' /W*), multiplying this singular term; as W* — > the singularity 
rapidly becomes weaker. Our approximation of the SIM to O (b) , given 
by (15), agrees very well with the SIM computed numerically by Lorenz 
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for b = 0.5. He finds Z* « 0.4PU* 3 for small W* for an initial condition 
corresponding to m — 0.81, while our equivalent result is Z* ps mbW* 3 = 
0.4051U* 3 . 

3.3 Divergent power series and periodic solutions on 



We now make a few observations about the expression (15) for the section 
through the SIM. Undoubtedly corrections will be made to Z* at orders b 3 , 
b 5 , . . . , but for our present purposes we assume that the essential features of 
the curve Z*(W*) are captured by (15) alone. 

Firstly, we recall the result that, in general, centre manifolds (and sub- 
centre manifolds, as are at issue here) are not unique. A given system may 
have an infinite number of centre manifolds, each differing from any other 
by exponentially small terms (see Carr [6]; and Roberts [25, §2.1]). Each 
has an identical power-series representation. Therefore the existence of more 
than one slow manifold (Ai, and the SIM, and infinitely many others) for the 
atmospheric model L86 should come as no surprise. The notable feature of 
the exponentially small difference between M. and the SIM is its multiplica- 
tion by a term with infinitely many algebraic singularities, cot(K/W*). This 
frustrating feature results from resonances between the fast and the slow 
waves that arise in attempting to construct a subcentre manifold, rather 
than a centre manifold (see Sijbrand [29]). Frequently, when one computes 
the power series for a centre manifold, it is divergent, even though a centre 
manifold contains no singularities similar to those of a subcentre manifold. 
We wish, then, to emphasize that the existence of singularities in the SIM 
and the divergent power series for any slow manifold are separate issues. 

Secondly, (15) indicates that the SIM and the slow manifold M. differ 
by terms that are exponentially small as W* — > 0. To see this we use the 
fact that we have chosen V* = in our search for periodic solutions of (1), 
and so according to the slow-manifold model, Z* = Z(U*, 0, W*). Using 
Roberts' [26] expansion of Z(U, V, W), (2), we find 



This agrees with the result (15) when the exponentially small term is ignored, 
since U* 2 = mW* 2 . Indeed, further investigation reveals agreement between 
W*F(W*) and Z(U* ,0,W*) at all orders in \(U*,W*)\. 

Thirdly, let us examine more closely the result of using the slow-manifold 
model (2-3) to approximate the SIM of (1). We begin by noting that the 



the SIM and on M 
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construction of Ai yields expressions for X(U,V, W) and Z{U,V, W) of the 
form 

X = UVf(U, V, W), Z = bWg(U,V,W). 

It follows that a zero-crossing of either U or V necessarily implies a zero- 
crossing of X. Further, the evolution of the slow variables U, V, W on Ai 
is governed by equations of the form 

U = -VW(1 -b 2 g) 

V = UW(l-b 2 g) (17) 
W = -UV, 

so that U 2 + V 2 is invariant for this model, just as it is for the original system 
(L86). Consequently, the dynamics of (17) are two-dimensional and so its 
solutions are in general periodic. Since, by definition, the SIM is composed 
of the periodic solutions of (1), it might therefore appear that Ai is the 
SIM of (1). However, we know from our preceding remark that Ai and the 
SIM of (1) are not the same — they differ by exponentially small terms. This 
discrepancy is resolved by concluding that a slow manifold Ai, constructed 
by the functional relationship y = h(x), cannot be exactly invariant un- 
der evolution according to (1) — there will always remain exponentially small 
discrepancies. 

The SIM was proposed because it can be computed directly without re- 
course to (divergent) power series. However, the power series for the SIM 
is the same as for Ai, and is divergent. We emphasize that the divergence 
of its power series does not mean that the SIM does not exist, nor does it 
mean that the SIM is ill-defined, or "fuzzy". The utility of the divergent 
series is illustrated by the excellent agreement between an approximation to 
Z*(W*) based on a small number of terms for Ai, and the numerical result 
for the behavior of Z*(W*) as W* — > on the SIM. For divergent series, 
this agreement is good for small W* and for a fixed number of terms. In 
the next section we shall illustrate some of the dynamical consequences of 
the differences between Ai and the SIM. For most solutions the differences 
are slight, and result in small quantitative changes. In a practical compu- 
tation, we are required to select one of the many possible slow manifolds. 
We consider two strong advantages of Ai over the SIM to be the possibility 
of computing appropriately the projection of initial conditions onto Ai, and 
the explicit algebraic form y = h(x) of Ai, compared with the construction 
of the SIM through an ensemble of numerical solutions. 
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4 Normal form for L86 

In Section 2 we described the construction of a slow manifold, Ai, which 
approximates a subset of the possible solutions of (1), and on which the fast 
"gravity" wave variables are slaved to the amplitudes of the slow "Rossby" 
waves. More general solutions of (1), with significant independent fast wave 
activity, involve oscillations centreed on Ai. In this section we rewrite (1) 
so that all solutions are captured, not just the slow solutions on Ai, but 
where the new governing equations describe the evolution of five new slowly- 
varying quantities. This is a normal form calculation (see Guckenheimer & 
Holmes [14]; Arnold [2]), which has been sketched for the forced, damped 
system of LK87 by Jacobs [17]. In that case the variables U, V, W, X, 
Z apparently may be written as convergent power series in the new slow 
variables: for the model L86, where there is no damping, we shall see that 
the equivalent series are divergent. Jacobs notes that in the new variables 
a slow manifold takes a particularly simple form. We shall also see from 
the normal form calculation that there is a limit on our ability to initialize 
data in such a way that the long-time dynamics of the physically dominant 
slow "Rossby" waves are unaffected by the balancing procedure. That is, 
there is a limit on the agreement between the unbalanced and the balanced 
systems — the two forecasts inevitably diverge slowly. 

We begin by noting that for the uncoupled system, (1) with b = 0, we 
can identify five slowly-varying quantities: U, V, W, R and 0, where X = 
RcosQ and Z = i?sinO. (In fact, when b = it follows that R = and 
= 1.) Our aim now is to find for the coupled system, (1) with b ^ 0, 
five equivalent slowly varying quantities. We do this by making successive 
nonlinear changes of variables: to leading order in the calculation the slowly 
varying quantities will be just U, V, W, R and 0. We shall perform algebraic 
manipulations on the system of equations, but there is no reduction in the 
dimension of the system as there was in computing a slow manifold Ai; 
we expect the normal form to capture all solutions of (1). This point is 
the essential difference between the simplifying tools of invariant manifold 
theory and normal form theory. The first aims to reduce the dimension of 
a nonlinear dynamical system, while the second transforms the system to a 
canonical form. 

To make the normal form transformation we seek to write (1) in terms of 
slowly varying variables u, v, w, r and 9. Linearly, these will be identified with 
the original variables U, V, W, R and 0, respectively. Thus the evolution 
of u, v and w will describe predominantly the dynamics of the slow waves, 
while r and 9 represent the amplitude and phase of the fast waves. We 
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shall expand the "physical" variables (U, V, W, X, Z) as power series in the 
new variables u = (u,v,w,x, z), where the two variables x = rcosO and 
z = r sin 9 are closely related to the original fast wave variables X and Z. 
This power-series expansion introduces exponentially small errors; the new 
evolutionary system for u will have solutions that differ from solutions of (1) 
by terms smaller than any power of |u| as \u\ — > 0. Further, if for a practical 
computation we truncate the transformation at a finite number of terms in 
the power series, then larger, algebraic errors are introduced. 

Recall that x = (U, V, W) and y = (X, Z) so that (1) may be written as 
the following pair of equations, 



x 

y 



f(x,y) 

By + g(x) 



(18) 
(19) 



where 



f(x,y) 



B 



(-VW + bVZ, UW 
1 





1 







bUZ^-UV), 
(0,bUV). 



We make a change of variables, 

x = X + F{x,ri), 



y = rj + G(x,rj), 



(20) 



where we identify x — ( u , v i w ) an d f] = (x, z). Here F and G are 0(|(%, ?7)| 2 ) 
as I (Xj r 7) I ~~ > 0) so that linearly x ~ x an d y ~ TJ- Substitution of (20) 
into (18-19) gives the following equations for the evolution of the new vari- 
ables: 



OF 

x + j^-x 

ox 

. dG . 
V + ^-X 
OX 



dF 
drj 
dG 
dr) 



f( X + F,ri + G) 

Bij + BG + gix + F). 



(21) 
(22) 



We now assume that the nonlinear terms and the new variables' evolution 
may be written as the power series 



F n, 

n=2 

oo 

X ~ E Rn(x,v), 



G ~ E Gn, 

n=2 

oo 

^7 ~ Brj + Yl Snix,*]), 



(23) 
(24) 



n=2 



n=2 



where F n , G n , R n , S n = 0(|(%, rf)\ n ). We shall solve (21-22) at successive 
orders in \(x, tj)\; at each order our aim is to choose the four quantities F n , 
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G n , Rn, S n as simply as possible. We shall see below what the term "simple" 
actually means. It turns out to be possible at each order to require that x, 
r and 9 are slowly varying, that is, that they are independent of the phase 
of the fast waves, 9. 

For use in the calculation below, we note that by the chain rule 

da da dx ^ da dz da ^ da da ^ 
dO dx d6 dz dO dx dz drj 

for any function a(u, v, w, x, z). As an example of some of the algebraic 
details, consider the equation governing the slow wave components R 2 and 
F 2 which is, from (21-24), 

dF 2 „ dF 2 

R * = f>--toj- B * = f>--dF> 

where 

f 2 = (—v(w — bz),u(w — bz), —uv). 

Note that f 2 depends on 9 through z. We now choose R 2 to be independent 
of 9 and such that secular growth in F 2 is avoided. This requires 

R 2 = (—vw, uw, —uv), 

in which case 

dF 

-^ = (bvz,-buz : 0). (25) 

Now we choose F 2 as simply as possible, that is, F 2 = (—bvx,bux,0). We 
can add to this expression for F 2 an arbitrary function of u, v, w and r, 
while still satisfying (25), but for simplicity we decide that ^-averages are to 
vanish. A consequence of this decision is that after averaging over the gravity 
waves u, v and w reduce to U, V and W respectively. 

For general n > 2, 

n-l n-1 xi? 

~ " d9 ' 



n 



where f n denotes terms of order \u\ n in the expansion of f(x,y) (and g 
is defined similarly from the expansion of g), and C n is a multi-nomial in 
x = r cos 9 and z = r sin 9 which involves quantities known from the previous 
calculation of R m , F m , S m , G m for rn — 2, . . . , n — 1. We choose R n to be 
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independent of 6, in fact to be C n , where the overbar denotes the mean with 
respect to 9, so that no secular terms arise in F n . Then we choose F n as 
simply as possible (that is, so that F n = 0). Note that R n is uniquely defined 
if it is to be independent of 9, and if no secular terms are to arise in F n . The 
term F n is unique only once we have specified F n . Altering the value of this 
average simply alters the relationship between the slow co-ordinates x an d 
the original variables by an amount of order \u\ n . 

Now we consider the equation for the fast wave components S n and G n , 
for n > 2, 



S n — BG n H — — — 

30 



n-l 



dG r 



9n 



m=2 



dx 



1 dGr, 

dr} 



(26) 



We note that D n can be written as a multi-nomial in r cos 9 and r sin 9 
which involves quantities calculated at previous orders. If we now compute 
5(26) + d(26)/d6 we obtain 



d 2 \ 




(27) 



where we have used the identity B 2 + 1 = 0. Our aim now is to make r and 
9 slowly-varying by requiring that they be independent of 9, that is, 



r 
9 



rp(x,r 2 ) 



But rr = xx + zz and r 2 9 = xz — zx, and so 



x z 
—z X 



r) 



X z 




X 




p 


—z X 




z 







Thus 



T] 



X 

z 



X 

z 



X 



p 



= M 



P 



(2? 



Jacobs [17, Appendix B] has termed (28) a "preferred choice" for the form 
of rj: at this point it is clear, though, that (28) is a necessary condition for 
r and 9 to be slowly-varying. 

Equation (28) implies that S n is of the form 5 n (x,r 2 )?7, where S n is a 
2x2 matrix. In (27) we choose G n to match all terms in the right-hand side, 
except the components of r cos 9 and r sin 9 in D n (for which 1 + d 2 / d9 2 = 0), 
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which must be removed by S n . So let us now consider these terms in (26). 
We know from (28) that S n is of the form 



M 



^n(x,r 2 ) 



It is also straightforward to confirm that if the components of G n propor- 
tional to x = r cos 9 and z = r sin 9 are 



G 



(i) 



7ix + 722; 
5\x + S 2 z 



where the 71, 72, Si and S 2 are functions of x an d r 2 , then 



-BG® + 



09 



x z 
—z X 



72 + Si 

$2 - 71 



N 



72 + Si 
$2 - 7i 



We note also that the terms in D n proportional to r cos 9 and r sin 9 may 
in all generality be written as Ma + Nb, where a, b are vector functions 
of x an d r 2 . This decomposition is unique, and therefore (p n , ip n ) = a and 
(72 + ^1,^2 — 7i) = b. Thus we have specified S n = a uniquely, although 
G n retains two undetermined coefficients. The two degrees of freedom left 
at our disposal correspond to redefining r and 9 by amounts of order \u\ n . 



4.1 Discussion of the normal form 

We have used the algebraic programming system reduce to implement the 
procedure described above to compute the normal form, which gives as the 
first few terms of the expansion 



u 


~ u — bvx + 


\b 2 u 


> 2 


-x 2 ) 


V 


~ v + bux + 




z 2 


-x 2 ) 


w 


~ w + bz{u 2 


-v 2 ) 






X 


~ x — buv + 


\b 2 x 


y 


-u 2 ) 


z 


~ z + bw(u 2 


-v 2 ) 







Despite the apparent differences, this normal form is equivalent to that of 
Jacobs [17], with his damping and forcing set to zero, because we have chosen 
our variables Xi V differently. 
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We now consider the structure of the normal form equations, and try to 
construct their SIM. A careful examination of the co-ordinate transforma- 
tion (20) at each order n, approximated above, reveals it be of the form 

U = u + uAi - bvA 2 

V = v + vA x + buA 2 

W = w + A 3 (30) 

X = x + X(u, v , w) + v4 4 

Z = z + Z(u, v , w) + A 5 , 

where (X(u, v, w), Z(u, v, w)) is the slow manifold M. computed in Section 2, 
and Aj = Aj(u,v,w,x, z) = O (r) as r — > 0. The evolution of the slow 
variables is given by expressions of the form 

u = —vw(l — Bi(u 2 ,v 2 ,w 2 ,r 2 )) 

v = uw(l — Bi{v 2 , u 2 , w 2 , r 2 )) 

w = -uv(l-r 2 B 2 (u 2 ,v 2 ,w 2 ,r 2 )) (31) 

r = ruvwC (u ,v ,w , r ) 

9 = 1 + D(u, v, w,r 2 ). 

The equations that govern the evolution of the four variables (u, v, w, r) 
representing amplitudes in (31) are independent of the fast gravity- wave 
phase variable, 9. We therefore examine first the four- dimensional sys- 
tem (31a-d) that results from ignoring the ^-equation, (31e). There are two 
invariants of the system (1), each independent of the gravity- wave phase: 
U 2 + V 2 and V 2 + W 2 + X 2 + Z 2 (L86). Similarly, these quantities are invari- 
ant under evolution of the normal form equations (31a-d), which therefore 
have two-dimensional dynamics. Solutions are therefore in general periodic. 
Re-introducing ^-variations, we see that solutions of the normal form (31) 
are in general quasiperiodic, with one frequency iO\ arising from (31a-d), and 
a second io 2 from (31e). Singly periodic solutions occur if io\ and uo 2 are ra- 
tionally related, or if r = 0; these solutions form the SIM of the normal form, 
which we construct below. The original system (1), however, not only has 
periodic and quasiperiodic solutions, but also has aperiodic solutions (LK87). 
Thus the normal form and (1), though quantitatively nearly identical, have 
qualitatively different dynamics. 

A further indication of the differences between solutions of (31) and (1) 
is the existence of some heteroclinic orbits in the former which are absent in 
the latter. To see this, we follow LK87 in considering the Hadley solutions 
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P V (F) of (1): namely the fixed points (U, V, W, X, Z) = (0, F, 0,0,0). (The 
same considerations will apply, with appropriate changes, to the solutions 
Pu(F): (U,V,W,X,Z) = (F, 0,0, 0,0).) Computations by LK87 indicate 
that gravity waves arise on the unstable manifold Ul{F) of each Hadley so- 
lution Py(F), except at isolated values of F. In these exceptional cases, 
Ul(F) is heteroclinic to Fy(F) and Py{— F), that is, Ul(F) is asymptotic 
to Pv{F) and Pv{—F) as t — > — oo and t — > oo, respectively. The near- 
heteroclinic behavior associated with almost all Hadley solutions causes so- 
lutions of (1) that pass through a neighborhood of the Hadley solutions to 
be aperiodic: solutions that remain away from the Hadley solutions are in 
general quasiperiodic, except on the SIM, where they are periodic. However, 
for the system (31a-d) the existence of two invariants and the symmetries 
of the system in any of the (u,v,w) co-ordinate planes constrain all Hadley 
solutions Pv(F) to be connected to their opposites Py(—F) by heteroclinic 
orbits. No aperiodic solutions of the normal form (31) exist. In summary, 
the dynamics of (31) may differ qualitatively from the dynamics of (1), par- 
ticularly for solutions that pass close by the Hadley solutions. The difference 
is due to the essentially two-dimensional nature of the normal form which 
occurs when 9 is decoupled, compared with the essentially three-dimensional 
nature of (1). The differences in the behavior of solutions of (1) and (31) 
that do not pass close to the Hadley solutions is small. 

A final illustration of the differences between solutions of (31) and those 
of (1) is given by the periodic solutions, which form the SIM. The SIM of (1) 
was described in Section 3 — to construct the SIM of (31) we first note that 
if we set r = at time t = then by (3 Id) r = for all time, so the plane 
r = is an invariant slow manifold of (31). Further, solutions with r = are 
in general periodic (because then the variations in 9 are irrelevant when the 
solution (U,V,W,X,Z) is reconstructed from (30)). Setting r = in (30), 
we see that the invariant manifold r = is precisely Ai, as calculated in 
Section 2. Some other periodic solutions of (31) occur for non-zero initial 
values of r — however, these have significant gravity-wave activity and they 
exist when the gravity waves happen to be exact harmonics of the Rossby 
waves. They give rise to manifolds of resonant solutions that intersect Ai, 
akin to the lines W* = K/ikix) for the uncoupled system (1) with b = 0. 
The slowest invariant manifold of the normal form (31) is therefore A4, and 
has no singularities. However, there are infinitely many "resonant" branches, 
which intersect A4 as shown in Figure 3. 

The relation between the normal form calculation of this section and the 
calculation of the slow manifold described in Section 2 is most easily seen if 
we assume that the new variables x are chosen, as we have done, by setting 
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F = 0. Then the transformation of the Rossby wave amplitudes is of the 
form 

x = X + F(X, V)=X + rF ix, V), 

so that by setting r = (that is, F = rj = 0) the new slow variables 
and the old become identical. To see the significance of setting r = we 
consider (21-22), which become 

x = f(x,h(x)) 

dh 

—x = Bh(x)+g(x), 

where h(x) = G(x, 0). These are precisely the equations that govern x and 
y = h(x) for an invariant slow manifold (Carr [6]). So by setting ry = 
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in (20) we recover Ai (implicitly, through y = h(x)). As a consequence, we 
expect the normal form sums in (23) and (24) to be divergent. 

There is a strong parallel between the method of normal forms and the 
method of averaging. (The latter method was suggested for the geostrophic 
approximation by Bill Dewar in private communication.) An alternative 
approach to the approximation of solutions of (1) is to assume that the gravity 
wave oscillations occur on a separate, faster time-scale than the evolution of 
the other variables, which are assumed to vary on the slow time-scale r. 
We then replace the operator d/dt in (1) by the operator djdr + d/d8, and 
write each dependent variable in (1) as a function of the new variables u{t), 
v(t), w(t), a(r) and the fast time 9. The notation for the new variables is 
as for the normal form calculation, except that the gravity waves now have 
a slowly- varying complex amplitude a(r). The two-time approach proceeds 
with similar algebraic steps as the normal form calculation, in particular at 
each stage an average over the fast time-scale 9 is taken in order to find terms 
in the slow evolution equations. 

4.2 Projection of initial conditions and the normal 
form transformation 

If the slow manifold had been a centre manifold, 2 characterized by an expo- 
nential approach of nearby trajectories to A4, then a rational principle can 
be used to compute the balanced initial data (Roberts [26]). The principle 
is that the long-term evolution of the original full system from its initial 
data must approach exponentially quickly the slow evolution on the centre 
manifold from the balanced initial data. This principle leads to an analysis 
whose framework is identical to that given in Subsection 2.1. However, the 
balanced initial data can be found immediately from the normal form of the 
equations, as we now demonstrate. Continuing the discussion for the case 
when A4 is a centre manifold, we find that a normal form transformation (20) 
can be chosen so that the new variables evolve according to equations of the 
form 



Thus r) = describes the centre manifold of the long-term evolution and, 
because the eigenvalues of B are negative, rj decays exponentially quickly to 

2 This is the case when damping of the gravity waves is incorporated into the model 
(Cox & Roberts [7]). 




(32) 
(33) 
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(as is the case when the fast gravity waves are damped). Thus (20) gives 
the following parametric description of the centre manifold in the original 
variables 

x = X + F(x,0), y = G( X ,0). 

However, from equation (32) we see that the evolution of the slow variables 
X is unaffected by the fast variables rj. Thus all solutions starting from 
points (x, v) with the same value of x have precisely the same x evolution, 
and because they all asymptote exponentially quickly to the centre manifold, 
so they all have the same long-term evolution. The hyper- surfaces in (x, y)- 
space that are described as rj is varied in (20) with fixed x form the isochronic 
manifolds introduced by Roberts [26] . Hence, the initial data (%, rj) are bal- 
anced to the centre manifold simply by setting rj = (as in a linear analysis). 
In terms of the original variables, the initial values (xq, Vq) are balanced by 
finding the corresponding normal form variables (xo, Wo) fro m (20), and then 
using 

x' = Xo + F(Xo,0), y'o = G(XoiO) 

as the appropriate initial conditions. From (x' Q , y' ) the evolution is slow and 
has precisely the same long-term dynamics as from (xq, y ). 

However, the "quasi- geostrophic" slow manifold M. is a subcentre man- 
ifold; it does not attract neighboring solutions exponentially, and the ar- 
guments above do not directly apply. The reason is that for a subcentre 
manifold the influence of the fast waves cannot be removed entirely from the 
evolution of the slow variables. This may be seen in (31) where the functions 
Bi and B 2 , which govern the evolution of the slow waves, also depend upon 
the amplitude of the fast waves, r: for a subcentre manifold (32) must be 
replaced, in general, by an equation of the form 

X = R(x,r 2 ). 

That is, the evolution of the slow variables cannot be entirely decoupled from 
the fast variables. Thus a simulation with both Rossby waves and gravity 
waves present need not be equivalent, over a long time, to any of the possible 
solutions with purely slow Rossby waves. 

This feature of subcentre manifolds is displayed in the simple system 

u = z 2 , x = —z, Z = X. 

These equations have the exact normal form, upon the transformation \ = 
u + xz/2, x = r cos 9 and z = r sin 9, 

X = r 2 /2, r = 0, 9 = 1. 
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Here the slow manifold is simply x = z = (r = 0), on which \ evolves 
according to x = 0. Thus the long-term evolution on Ai is trivial: x is 
constant. However, if there are any fast waves present, r ^ 0, then \ drifts 
according to % = r 2 /2, but there is no corresponding solution for \ on A4. 
The same is true for the normal form (31) of the model L86 — at fourth order 
in |x| we cannot avoid introducing gravity- wave terms proportional to r 2 
into the evolution equations for the slow Rossby waves. Therefore solutions 
on and off the slow manifold M. cannot share the same values of u, v and w 
for all time; there are inevitable discrepancies which become significant on a 
timescale of O (r -2 ). 

It is now clear that the balancing procedure, the projection of initial 
conditions described in Section 2, which is linear in the fast gravity wave 
amplitude r, cannot be improved to be correct to O (r 2 ). 

5 Conclusions 

A significant part of our understanding of atmospheric dynamics rests on 
the concepts of quasi-geostrophy (Gill [12, Chapt. 7]). The fundamental 
concept is the separation of the dynamics into waves of two time-scales: the 
slowly evolving Rossby waves and the quickly evolving gravity waves. We 
have described how a slow manifold, formally composed of the ensemble of 
slow Rossby waves, may formally be written in the form y = h(x), with h 
developed as a power series in the slow variables x. This is similar to the 
scheme proposed by Baer & Tribbia [3], and is equivalent to that of Vautard 

6 Legras [32]. As L86 demonstrates, the power series for h is divergent, so 
that successive approximation schemes at first appear to converge to a slow 
manifold, but after the inclusion of sufficiently many terms they begin to 
diverge. Nevertheless, the divergence of its power series does not indicate 
the non-existence of a slow manifold. Each successive approximation to M. 
captures the dynamics of L86 to within a higher power of x, but in a smaller 
radius around the origin. Unfortunately for our purposes, the dynamics on 
M. (determined exactly, but not from its asymptotic series) differ from those 
of L86 by an amount that is smaller than any power of X clS X ^ 0. The 
dynamics of the closed set of evolution equations (3) for the slow variables 
on M. are genuinely slow: no gravity waves develop. But M. is not quite 
invariant, so the apparent slow dynamics on Ai do not reflect genuine slow 
dynamics in L86. In particular M. appears to be filled with periodic solutions, 
whereas the full system has a manifold of periodic solutions (Lorenz' sim) 
which is different from A4 and which also contains singularities. 
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The issue of how to balance given data measurements by projecting them 
onto a slow manifold is a complicated one for realistic situations — the most 
appropriate projection may depend on the relative quality of the various mea- 
surements, and on the influence of different measurements on the projected 
initial point (Daley [9, 10]). We have described a method of projection that 
is certainly well-suited to models where the initial data are known exactly, 
such as L86. The criterion we have proposed for selecting an initial point 
on M. seems appropriate for numerical weather prediction: namely, that the 
behavior of the balanced system should correspond for all time to that of 
the unbalanced system (just without the fast gravity-wave activity). This 
ensures that the essential features of the forecast are unchanged by the ini- 
tialization process. However, we have also shown that due to resonances, 
as exhibited by the normal form transformation (Arnold [2]), the initializa- 
tion can be carried out only to first order in the fast wave amplitude 3 . We 
have illustrated the practical utility of our proposed initialization procedure 
with numerical integrations (Figure 1) that show the superior accuracy of 
the forecast from appropriately initialized data. 

Certain difficulties associated with computing a slow manifold have led 
researchers to question the existence of such a manifold, or to label the con- 
cept "fuzzy". A slow manifold for L86 is a subcentre manifold (Kelley [18]; 
Sijbrand [29]), and is not unique: just as for all low- dimensional models 
(Roberts [27]) there are infinitely many slow manifolds in which the fast vari- 
ables are specified in terms of the slow variables. It is this non-uniqueness 
that accounts for the "fuzziness" in the concept of the slow manifold: the 
partial differential equation that governs a slow manifold for a given problem 
is perfectly well-defined, as is each manifold. It is just that the construction 
of low- dimensional models has infinitely many solutions, differing by expo- 
nentially small terms. It is possible to select one distinguished slow manifold 
by assuming that the fast variables may be written as power series in the 
slow variables. However, often, as is the case for L86, the power series for 
the slow manifold proves to be divergent, so that the series is asymptotic. 
As an alternative, the numerical construction of a slow manifold (Lorenz' 
Sim) was conceived to avoid the divergence problem. The SIM is defined as 
containing all the periodic solutions. However, in the light of our previous 
comment that slow manifolds differ by exponentially small terms, the SIM 
must share the divergent power series of all slow manifolds for L86. This 

3 The resonances that force certain terms to occur in the evolution equations for the 
normal form variables are a consequence of the linear dynamics (Arnold [2]), and are quite 
distinct to the resonances between the periods of the (fully nonlinear) Rossby and gravity 
waves that induce singularities in the slowest invariant manifold. 
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does not imply that the SIM does not exist, nor does it suggest that the SIM 
is somehow "fuzzy". The SIM exists; it just has a divergent power series. The 
asymptotic series of all slow manifolds for a given system are identical — the 
slow manifolds differ by sub-dominant terms (Bender and Orszag [4]). 

We contrast our point of view with that presented by Warn & Menard [33] , 
relating to the slow manifold for the model L80. They argue that a slow man- 
ifold must attract neighboring solutions in order to be useful in applications. 
Since almost every numerical solution of L80 involves some level of fast "grav- 
ity" wave activity, they conclude that an attracting invariant slow manifold 
does not exist. However, it is important to recognize that a slow manifold 
need not be attracting in order to serve as a useful centre for the dynamics 
of the system. The subcentre manifold we have discussed for L86 does not 
attract neighboring solutions, yet it has dynamics that approximately corre- 
spond to the full initial-value problem, except that the fast gravity waves are 
absent. Warn and Menard argue that the slow manifold should be replaced 
by a more general "fuzzy" balanced set because their scheme to compute a 
slow manifold encounters asymptotic series, and because in their approxima- 
tions to the slow manifold, gravity wave activity is observed. We have argued 
that a slow manifold is necessarily "fuzzy" because of its exponential non- 
uniqueness, which is inherent to invariant manifolds of the centre-manifold 
family. Each individual invariant manifold is well-defined. The "fuzziness" 
of the set of manifolds is an entirely separate issue from that of the diver- 
gence of the power series for each slow manifold. As Lorenz has illustrated 
for the SIM, one may compute a slow manifold numerically, even when it has 
a divergent series. 

We cite two other examples of common and useful approximations that 
correspond to the use of a subcentre manifold, even though such a mani- 
fold does not attract neighboring solutions. The first is the incompressible 
approximation in fluid mechanics, where the fast variables represent sound 
waves and are neglected to leave the slow evolution of incompressible flow. 
The second is beam theory in elasticity, where rapid internal elastic vibra- 
tions are the fast variables and where large scale deformations are the slow 
variables (Roberts [28]). The fast "ringing" modes are eliminated in the tra- 
ditional approximations of beam theory. The concept of a subcentre manifold 
enables one rationally to extend these traditional approximations to higher 
order, to derive appropriate initial and boundary conditions for the model, 
and to treat forcing appropriately. 
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A Contour integration 

In order to compute Z%, we wish to evaluate the integral 

I = U* 2 [ T ° sin(T - r) cn(W*r) sn(W*r) dr, 
Jo 

which appears in (14). We first expand the trigonometric function, then 
change variables by setting t = W*t, so that 

I = mW* ( I c sin T - I s cos T ) , (34) 

where 

r K t 
h = J Q cos — sn(t) cn(t) dt, 



and 

rK 



f n t 

I s — / sin sn(t) cn(t) dt. 

Jo W* 

To evaluate I s , we consider the contour integral 

I r = [ S(t; W*) dt, (35) 



r 

where 

S(t;W*) =e it/w *sn(t) cn(t), 

and where T is the rectangle with vertices at —K, K, K + 2iK', —K + 2iK' as 
shown in Figure 4. We denote the straight line segments joining these vertices 
from — K counterclockwise by Ti, T2, T^, T^, respectively, with corresponding 
integrals I±, I2, h, h- Then 

-K 



h = J e it/w *sn(t) cn(t) dt = 2il s . 
The integral along the parallel side of T is 

J 3 = J~ K e i{2iK ' +t)/w 'sn(2iK' + t) cn(2iK' + t) dt 

-IK' /W* t 

— e ±\. 

The other integrals along the vertical sides of the rectangle give 

r 2K' 

h = / e iiK+it)/w * sn( K + it) cn(K + it) idt 
Jo 

r2K' 

= m\ /2 e iK/w * / e~ t/w *nd(t\ mi ) sd( t\ mi )dt, 
Jo 
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-K + 2iK' 



K + 2iK' 




Figure 4: Paths of integration in the complex plane. 



where mi = 1 — m, and 



.T* 
1 2- 



The contour Y encloses a double pole of S(t; W*) at t = iK, and so 

I r = I 2 - I* + 2i(l + e- 2 ^' /w/ *)/ s = 27rzResidue(S(t; = itf')- 

To compute the residue we expand S(t; W*) near the pole to give 

t — iK'\ ( m" 1 / 2 \ ( —imr 1 ! 2 ' 



S(t; W* 



-K'/W* 1 + 



W* 



t-iK' \ t-iK 



T + 0(1) 



where the singular behavior of the Jacobian elliptic functions is (Table 16.7 
of Abramowitz and Stegun [1]) 



sn(t) 



m 



-1/2 



t-iK 



7 , cn(t) 



—im l l 2 . 
^ as t — > iK . 



t-iK' 



Thus 



Therefore 



Residue(S(t; W*);t = iK') 



1 



mW* 



-K'/W* 



I * = rr^Sech 



m 



1/2 



2mW* 



W* 1 + e ~ 2K '/ w * Jo 



2K' 



-t/W* 



nd(t|mi) sd(t|mi) dt sinT 
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Through successive integrations by parts we find that the integral in this 
expression has the asymptotic series 



r 2K' fi- 



r2K' 



1 d 



1 oo 

(l + e - 2K '/ w ') ^* 2n cd (2n) (0). 

m n=l 

Thus we may write 

I s — sech^- + ^- V W- 2n cd^(0) sin T . 



n=l 



Now we note that the integral I c has the asymptotic series given by 

rK I 

I c = cos — sn(t) cn(t) dt 

f K t ( 1 d , .A , 

= / cos — — rdn(r) ore 

Jo W* \ mdt y 7 

1 oo 

~ E W* 2n (-l) n " 1 (dn (2n) (K) cosT - dn (2n) (0) 



71=1 



and that 



dn (2n) (^) = mi /2 nd (2n) (0) = (-l) n mi /2 cd (2n) (0) 



because nd(it|m) = cd(t\mx). Therefore, by substituting the expressions we 
have calculated for I s and I c into (34), we find 



mW* 



1/2 



m 
m 



J-/^ OO 1 oo 

i_ ty* 2 - c d (2n) (0) cos T Y, ^* 2n dc (2re) (0) 

12 71=1 m 71=1 



sinTn 



7T 



sech-^ + ^- £ Vr 2n cd< 2n) (0) sin T 



2mW* W* m 



71=1 



cosT 







1 oo 

-]T vr 2n dc (2n) (o) 



71=1 



iT' 

sin T — — sech— — cos T 



This is the expression we desire for /. The sum that appears is divergent, 
because the Taylor expansion 



OO 1 

dc(W*) = J2 7^-T7^* 2n dc (2n) (0) 
n=0 \^ n )- 
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has a finite radius of convergence (equal to K). It follows now from (14) that 
Zl = —I I sinT 

= J2 ^* 2n+1 dc^ (0) + ^sech^ cot T . 

n=l 1 W 
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